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' In two papers we proposed a continuum model for the dynamics of systems of self propelling 

particles with kinematic constraints on the velocities and discussed some of its properties. The 

^ , model aims to be analogous to a discrete algorithm used in works by T. Vicsek et al. In this paper 

O ' we derive the continuous hydrodynamic model from the discrete description. The similarities and 

' differences between the resulting model and the hydrodynamic model postulated in our previous 

, papers are discussed. The results clarify the assumptions used to obtain a continuous description. 
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I. INTRODUCTION 



Systems of self propelling particles (SPP) are widely represented in the nature. Flocks of birds, schools of fishes, 
swarms of ants, groups of bacteria etc. are examples of such systems [H 0; S I3- The observed complex behavior of 
^ r these systems is far from completely understood. That is why this phenomenon is of great interest. Pattern formation 
I in the systems of SPP is caused either by the presence of external fields (e.g. temperature gradient, food density 
c/j gradient, chemotaxis phenomenon etc.) or by kinematic constraints which are imposed on the motion of the particles. 
O In this second case the clustering of SPP is driven by internal dynamics of a nonpotential character. The tendency 
c/3 \ of the particles to align their velocities with their neighbors is a crucial element in the mechanism of the emergence 
J>~^' of the coherent motion of the SPP. Numerous attempts have been made to find a model describing the collective 
behavior of self propelling particles. One may distinguish two main directions of research: the numerical simulations 
I (discrete description) and the hydrodynamic (continuous) approaches. 

The first numerical model (discrete algorithm) for coherent motion of SPP was proposed by T. Vicsek et al. [H]. 
T-H In their work a simple kinematic updating rule for the directions of the velocities of the particles was proposed 
^ and numerical evidence was given for a transition from a disordered state to ordered collective motion at low noise 
' amplitude and high density values. They seem to have been the first to realize that flocks fall into a broad category 
of nonequilibrium dynamical systems with many degrees of freedom and noted the analogy between flocking and 
I ferromagnetism. The model of T. Vicsek has become the "minimal" simulation model in the "modern era" of the 
^-H ' study of flocking. Extensions of T. Vicsek's model have been proposed considering particles with varying speeds 
\^ ' and different types of noise, by inclusion of external force fields and/or intcrparticle attractive and repulsive forces 

o: saii; _ 

(/2 ' Properties of T. Vicsek's model were also investigated from a mathematical point of view. In [lO| the spontaneous 
O , emergence of ordered motion has been studied in terms of so called control laws using graph theory. Generalizations 
' c/j ' of the control laws were considered in [ll|, ■ In it was shown that the organized motion of SPP with the control 
J>^| laws depending on the relative orientations of the velocities and relative spacing, can be of two types only: parallel 
■ and circular motion. The stability properties of these discrete updating rules (including the T. Vicsek's model) and 
the dynamics they describe were considered using Lyapunov's theory in [Tol. [Til. [Tsl. [T^ . 

The work by T. Vicsek et al. also gave an impulse to the development of continuous approaches. Here one may 
distinguish two classes of approaches. 
^ . The first class consists of the models which are usually based on the analogy with the Navier-Stokes equation, 
' where terms describing the self propelling nature of the system are added. In [1^ the pressure and viscous terms are 
incorporated into the model side by side with the dri ving force and the friction caused by the interaction with the 
environment. The inclusion of additional terms in [TgI. Il7j| was done based on symmetry consideration. The attempt 
to derive such phenomenological equations from the kinetic equation was made recently in p^ . 

The second class contains models which describe the swarming behavior of SPP by inclusion of attractive and 
repulsive interactions. The model, based on the diffusion- advection-reaction equation with nonlocal attractive and 
repulsive terms, is suggested in [T9| in order to describe the swarming phenomenon. Their model gives a compact (with 
sharp edges) aggregation of SPP with a constant density as a result, which according to the authors is biologically 
reasonable. 

Another continuous model for the behavior of the living organisms with nonlocal social interactions is proposed in 
[20| . There the kinematic rule for the velocity field contains the density dependent drift and the nonlocal attraction 
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and repulsion contributions. For a 2-dimensional case of incompressible fluid with the motion of the particles being 
normal to the population gradients, the flow of a constant density with a compact support is obtained. 

In two papers [2lL |22| we proposed a hydrodynamic model, which can be considered to be the continuum analogue 
of the discrete dynamic automaton proposed by T. Vicsek et al. 0], which we further will call the CVA model or 
algorithm. We constructed our model on the basis of the physical properties of the discrete CVA model, namely the 
conservation of the number of particles and the kinetic energy. The discrete configuration updating rule used by T. 
Vicsek et al. changes only the direction of the particle velocities but keeps their absolute value constant. In their 
algorithm the number of particles is constant as well. 

In this article we obtain the continuous description by coarse graining the discrete CVA algorithm. In this respect 
our present paper is meant to be a link between two existing groups of approaches: discrete and continuous. The 
importance of this analysis is that it clarifies which of the continuous models we proposed is closest to be the continuum 
analog of the CVA model. 

In Section 2 we will start with a rule for the velocities formulated by T. Vicsek et al. and obtain a discrete equation 
of motion for each particle. We introduce angular velocities associated with the rate of change of the direction of the 
linear velocity of the particles. These angular velocities contain the information about the nonpotential interactions 
between a given particle and its local surrounding. We derive an expression for the angular velocities in the continuous 
time description. We show that to a first order in the velocity difference between the steps the angular velocity for 
particle i depends on the average velocity in the neighborhood of the ith particle and its rate of change. 

In Section 3 we obtain the continuous description, with a conserved kinetic energy and number of particles, using a 
coarse-graining procedure. We obtain the angular velocity field that follows from the CVA 2-dimensional model and 
compare it with the angular velocity fields we proposed in our first paper. It turns out that there are similarities and 
differences. Both the continuous description that follows from the CVA model and our continuous model are and have 
been shown to give stationary linear and vortical flow fields. The description of such flow fields is one of the aims of 
the model. A discussion is given and concluding remarks are made in the last section. 



II. CONTINUUM TIME LIMIT 



In this section we derive the equation of motion in continuous time from the CVA algorithm. In their work ^ the 
collective behavior of self propelling particles with respect to a change of the density and the strength of the noise 
was investigated. In our analysis the noise will not be considered. We focus on the systematic contribution. In our 
first paper we discussed how noise can be added in our approach. 

The ordered motion of self propelling particles in jsj is described by the CVA rule, according to which at each 
discrete time step (labeled by n) the ith particle adjusts the direction of its velocity (n) to the direction of the 
average velocity (n) in its neighborhood. The average is calculated over a region with a radius R around a given 
particle. Using this radius we will call particle densities small compared to R^'^, where d is the dimensionality, small. 
When the density is larger we call it large. The CVA rule implies that 



Vi {n + 1) X Ui (n) — 0, yi,n, 
where the absolute value of the velocity of each particle is assumed to be constant , i.e. 

\vi{n + 1) 1 = 1 (n) 1= Vi . 

Together with Eq. ([T|) it follows that 

Vi (n + 1) = Vi (n) , where | (n) |= 1 . 



(1) 
(2) 
(3) 



Using the fact that (ra + 1) — (n) is perpendicular to v.^ (n + 1) + (n), given the validity of Eq. it can be 
shown that 



Vi (n+l) ~ (n) = [vi (n) x (n + 1)] x 



V, (n -I- 1) -I- (n) 
1 + Vi (n) ■Vi{n+ 1) 



(4) 



where (n) = v.; (n) / Vi is a unit vector in the direction of the velocity (n) . 

It is important to realize that there is a difference between low density regions and high density regions. In 
high density regions the velocity of the particles is updated at every step. In the low density regions the average 
of the velocity of particles around and including particle i is equal to the velocity of particle i. It follows that 
Ui [n) = Wi {n) /vi. As a consequence (n + l) = (n). Substitution in Eq. ([4]) gives the equality zero equal to 
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zero. The important conclusion is that in the low density regions the particles do not change their velocity. We will 
come back to this point when this is relevant. 

In order to obtain a continuous description as a function of time, we assume the steps to be small so that 

9, + (n) |<Cl. (5) 

One may then write Eq. ^ to first order in the velocity difference as 

Vi (n + 1) - Vi (n) = [vi (n) x v,; {n + I)] x = [v^ (n) x (n)] x (n) . (6) 

As we are interested in the rate of change of the velocity we divide this equation by the time step duration r. This 
gives 

Vi (n) X Ui (n)' 



(n + 1) ~ Vj (n) 

T 

where 



X Vi [n] = LJ^. X V, , (7) 



1 ^ 

- X Ui (8) 

T 



is an angular velocity associated with the velocity vector v.; . 

In view of Eq. ([5]) the left hand side of Eq. ([7]) gives the continuous time derivative. In other words we may introduce 
the following definition: 

V, in) ^ -^ir^+^)~-dn) 
r 

Using that (n) = (n + 1) it follows from Eq. fT]) that 

Ui (n + 1) - Ui (n) = TWu, (n) x u^ (n) , (10) 
where the angular velocity u:^. (n) corresponding to the average velocity u^ (n) is defined as 

^^ui (n) = u:^^ (n + I) . (11) 

It can be shown that 

t^u. {n) = - [vj (n + l) X Ui {n + 1)] = - [u^ (n) x u^ (n + 1)] = u^ (n) x u^ (n) (12) 

T T 

where 

. , ^ u, (n + l) - u, (n) 

Ui (n) = . (13) 

T 

Furthermore one may show that 

LJu, (n) - Wv, ("-) = (n + 1) X Vi (n) = r u^ (n) x Vi (n) , (14) 
where the second order derivative is defined by 

(n) = \ W^ {n + 2)- 2vi (n + 1) + Vi (n) 1 . (15) 

Combining Eqs. (fT2|) and (fT4|) results in 

a^v, {n) = uj^^ (n) - rUi (n) x v., (n) = u, (n) x li^ (n) -tu, (n) x Vi (n) , (16) 

which to first order in the velocity difference implies that for the angular velocity associated with the particle velocity 
Vi we obtain the following expression: 

(n) — u)ui (n) — TUi (n) x u^ (n) = Ui (n) x li^ (n) . (17) 

The second equality follows from the fact that the second derivative iii (n) is parallel to Ui (n) to first order in the 
difference. 

Replacing n by the time t the resulting equation of motion becomes: 

= LJ^^{t) X v,(t) = [u,{t) X U,(t)] X v,(t). (18) 

This equation is continuous in time and is derived from the discrete CVA rule using Eq. ([5]). In order to obtain 
equations for the velocity and the density fields, which are continuous in space, we will coarse-grain Eqs. p7p and 
in the next section. We note that both c?Vi (t) /dt and iii (t) are zero in the low density regions. 
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III. CONTINUOUS TRANSPORT EQUATIONS 



In this section we introduce the averaging procedure for the discrete model. We derive continuous equations for 
the vefocity and the density fields by averaging the discrete equations. 

In our continuum model the number density, n(r,t), satisfies the continuity equation, 



dn{r, t) 
dt 



+ div (n(r,i) v(r,<)) = 0. 



(19) 



The dynamics of the velocity field is such that the kinetic energy density is conserved, d\\-(r,t)\'^ /dt — 0. The time 
derivative of the velocity field is therefore given by 



(20) 



where ijJ (r, t) is some angular velocity field. We will discuss how to obtain this angular velocity field in terms of the 
velocity and density fields below. 

In the CVA algorithm the direction of the average velocity in the neighborhood of particle i is given for the 
continuous time description by 



{t) = Y.H in, (t)) V, {t)\j2H (r„- (t)) V, (t) 



(21) 



where = jr^ — |. The dynamics of individual particles therefore reduces the difference between the direction of its 
velocity and that of the average velocity of the surrounding particles. H (r) is an averaging kernel, which we assume 
to be normalized, 



J H{r)dr = l 



(22) 



It has the characteristic averaging scale R, beyond which the kernel goes to zero Usually one uses for H a 

normalized Heaviside step function. 

In order to obtain a continuous description we define the average particle density (per unit of volume) and velocity 
fields by 



n(r,t)=^i/(r-r,(t)), 

j 

,(r,^)v(r,^) = ^i^(r-r,(^))v,(^). 



(23) 



Using Eq.(I221) in Eq.dni) for r = r^, it follows that (<) = v (r, (t) , t), where v (r^ (t) ,t)^v (r, (<) , i) / |v (r^ (t) ,t)\. 
Eq. can therefore be written as 



dt 



v{r^{t),t) X 



dv(r,(i),t) 



dt 



V. {t) . 



By evaluating the time derivative between the square brackets on the right hand side one obtains 

dvi{t) 



dt 



[9 (r, t) X (v, it) ■ Vv (r, t))^^^^^^^ x v,(t)+ 
dv (r, t) 



V (r, t) X 



dt 



X Vi{t) 



r=ri(t) 



(24) 



V (r, t) X [(v, (t) - v(r, t)) ■ Vv (r, t))]r=r,(t) x ^.(0+ 
dv(r,i)" 



V (r, t) X 



dt 



X Vi{t) 



(25) 
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In view of the fact that the gradient of the direction of the average velocity is of the first order and that the velocity 
difference is also of the first order, the first contribution on the right hand side is of the second order and can be 
neglected. Eq. ([25]) therefore reduces to 



dt 



v(r,t) X 



dv (r, t) 
dt 



(26) 



r=ri(t) 



Note that both dvi {t) /dt and dv (r^ (t) ,t) /dt are zero in the low density regions. When we now average this equation 
we can use the fact that the expression between the square brackets only depends on the coarse grained functions and 
therefore varies slowly over the range of the averaging function. 

Averaging the left hand side of Eq. (j26|) and using the continuity equation we obtain 

^ 5^ V, (t) i/ (r - r, (t)) - E (^) 4r ■ (^) ^i"^' (^)) = 



dt 



d {n{r,t)v (r,t)) d v-^ , , , 

+ ■ E ^« (*) H{v-v, (t)) 

i 

n (r, t) ^-^^^ -V-[n (r, t) v (r, i) ® v (r, t)] + 



V • ^ V,; {t) V, {t) H{v-Y, (t j) = n (r, t) 



dv (r, t) 
Jt 



V • 5] (v (r, t) - V, {t)) (v (r, t) ~ v, (t)) H {v - (t)) 

i 

dv (r, t) 



n{r,t) 



dt 



(27) 



where we neglected the term of the second order in the velocity difference. This term would give a small contribution 
to the pressure tensor. It is therefore also a contribution which is in the formulation of the problem assumed to be 
cancelled by the self propelling force. 

In the right hand side of Eq. we have 



E 



■ir,t) X 



dw{r,t) 



V (r, t) X 



dt 
dv (r, t) 



XV, {t)H{Y-r, (t)) = 



dt 



r=r>(t) 

X n (r, t) V (r, t) . 



This implies that the averaged equation of motion can be written as 



dv (r,<) 
di 



V (r, t) X 



dv (r,<) 
Jt 



X v(r,t) 



This gives 



dv (r, t) 



dt 



u: (r, i) X V (r, t) 



where to first order in the velocity difference 



u){Y,t) = v(r,t) X 



dv (r, t) 
dt 



v2 (r,t) 



V (r, t) X 



dv (r, t) 
dt 



(28) 



(29) 



(30) 



(31) 



where v (r, t) = | v (r, t) \ . 

Here we restrict our discussion by considering the 2-dimensional case in order to make a comparison with the results 
obtained in our previous papers. By evaluating the time derivative in Eq. pip we may rewrite this expression as 
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follows 



cj(r,t) 



v2 {r,t) 



v{r,t) X +(v(r,i).V)v(r,t) 



v2 (r,t) 



v2 (r,i) 



V (r, t) X [v (r, i) x rot v (r, t)] 



V (r, t) dv (r, t) 



v2 (r,t) 



at 



v(r,t) „v2(r,t) , , 



(32) 



This is the angular velocity field obtained from the discrete algorithm used by T. Vicsek et al. 

One may see that the continuous equation of motion, Eq. ((30|) . with the angular velocity derived from the CVA 
rule, Eq. ([5^ . can be written as follows: 



dv 



(1 - vv) 



5v 

'dt 



(1 



vv) • V— + (rotv) X V. 



(33) 



where 1 is the unit tensor. All three terms on the right hand side contribute to the co-moving derivative of the velocity 
which is orthogonal to the velocity field. 

In the low density regions one obtains, as has been pointed out a number of times, /dt — 0. As this is not so 
clearly visible in Eq. p3)) it is appropriate to replace this equation by 



dw 
'dt 



(1 - vv) • ^ + (1 - vv) • Vy + (rotv) X v 



(34) 



where f{n) is the density dependent factor, which arises due to coarse-graining procedure. In low density limit it is 
natural that f{n) ^ as n ^ 0. A more thorough analysis of this is beyond our present aim, however. 

Before comparing this expression to the one we used in Refs. [2l|, we first verify that stationary linear and the 
vortical solutions are solutions of Eq. (|33|) . In view of their stationarity the first contribution in Eq. (l33l) is equal to 
zero. For a linear flow v = vq e^; the other two terms and the left hand side of Eq. ([55)) are also zero. Stationary 
linear flow is therefore a solution. In case of stationary vortical flow, v = y^pir) e^p {(fi), the v • Vv term on the left 



hand side of Eq. ([55]) cancels the terms due to the second and the third term. The continuity equation, Eq. p^ . is 
satisfied for each density distribution which varies only in directions normal to the fiow direction. It follows that the 
continuous analog of the CVA model has stationary linear and vortical solutions. 

In our first paper [2l| we used an angular velocity field which was a linear combination of n (r, t) rot v (r, t) and 
Vn (r, <) X V (r, t). The resulting equation of motion becomes 



dv , 

— = Sin (rotv) X V -(- S2 (Vn x v) x v. 



(35) 



The first term is analogous to the third term on the right hand side of Eq. ([33| . Similar to the CVA model this 
choice leads to stationary linear and vortical solutions. The linear dependence of our choice on the density leads to a 
dependence of the stationary velocity field on the density distribution. We refer to [HI, [23] for a detailed discussion 
of these solutions. For a small density the right hand side of Eq. ([55)1 makes dv/dt negligible. This is similar to the 
behavior in Eq. ([34]). 

When one modifies the updating rule in the CVA model, as done in Refs. [3, [1, [^, this leads to a modification of 
the u) given in Eq. ([5^ . Similarly, the choice of uj we used in [2|[[23| can be modified to include such contributions. 
The freedom in the choice of uj in the continuous version of the CVA model is one of its strength. 



IV. CONCLUSIONS 

In our first paper [23| we constructed a continuous self propelling particle model with particle number and kinetic 
energy conservation. In this paper we addressed the problem to derive the continuous description from the discrete 
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model proposed by T. Vicsek et al. [5|. We were able to derive expressions for the angular velocity field used in the 
continuous model from the updating procedures used in their model. By coarse graining the discrete equations in the 
original model we obtained the angular velocity field used to give the co-moving time derivative of the velocity field 
in the continuous description. Modification of the updating rules in this model, as done in Refs. @,02@, results in 
modifications of the resulting angular velocity field. The angular velocity field used in our work [2ll.l22l| is one of such 
choices. One of the contributions in the continuous version of the CVA model is very similar to one of the contributions 
which we have postulated in our hydrodynamic model \21]. Both the continuous CVA model and our model lead 
to the linear and vortical flows of the self propelling particles observed in nature and obtained in simulations and 
continuum approximations. This shows that they are appropriate for the description of flocking behavior, which is 
one of the aims of the model. An interesting alternative choice of a; (r, t) in the continuous description is for instance 
Wca X V where ca is the concentration of an attractant. As was shown in A. Czirok et al. [7] this choice can be used 
to describe rotorchemotaxis . Note that the term Vn (r,i) x v (r,t), which we considered in our continuum model, is 
similar to the one considered by A. Czirok et al. when the concentration field is proportional to the concentration of 
the attractant. 

Our analysis shows that one may coarse grain the discrete updating rules and obtain the corresponding continuous 
description. This makes a direct comparison between discrete and continuous descriptions possible. For our own work 
it was found that our continuous description was similar to the continuous version of the original CVA model but not 
identical. The analysis in this paper makes it possible to extend our work on the continuous description such that it 
is either closer or more different from the original CVA model. 
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